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The Olami-Feder-Christensen earthquake model is often considered the prototype dissipative 
self-organized critical model. It is shown that the size distribution of events in this model results 
from a complex interplay of several different phenomena, including limited floating-point precision. 
Parallels beween the dynamics of synchronized regions and those of a system with periodic boundary 
conditions are pointed out, and the asymptotic avaanche size distribution is conjectured to be 
dominated by avalanches of size one, with the weight of larger avalanches converging towards zero 
as the system size increases. 
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Self-organized critical (SOC) systems |l], || are ex- 
tended systems that receive a slow energy input and have 
fast dissipation events ("avalanches" or "earthquakes"), 
the size distribution of which is a power law. Such a scale 
invariance can result quite naturally in systems that sat- 
isfy a local conservation law however, the mechanisms 
leading to SOC in non-conserved systems are not yet 
well understood. The Olami-Feder-Christensen earth- 
quake model jl] is probably the most studied nonconser- 
vative SOC model. Nevertheless, the nature of its critical 
behavior is still not clear. Despite the simplicity of its 
dynamical rules, the model shows a variety of interest- 
ing features that are unknown in equilibrium physics and 
appear to be crucial for generating the critical behavior. 
Among these features are a marginal synchronization of 
neighboring sites driven by the open boundary conditions 
H , and the violation of finite-size scaling ](| f?J together 
with a qualitative difference between system-wide earth- 
quakes and smaller earthquakes Also, small changes 
in the model rules (like replacing open boundary condi- 
tions with periodic boundary conditions || , or introduc- 
ing frozen noise |l(|), destroy the SOC behavior. 

The model is a discretized and simplified version of 
the Burridge-Knopoff model of earthquakes and is 
defined by the following rules: At each site of a square 
lattice, a continuous variable z is defined that represents 
a local force. The force at all sites increases uniformly 
at constant rate. When the force z(i,j) at a site (i, j) 
exceeds the threshold value z c , which can be chosen to 
be z c = 1 without loss of generality, the force at this site 
is reset to zero, while all four nearest neighbors receive a 
force increment of az(i,j). The parameter a is the only 
parameter of the model, and it has a value in the interval 
(0,0.25). If a neighbor is lifted above the threshold, the 
force on its neighbors is increased according to the same 
rule, etc., until the "earthquake" is finished. The "size" 
of an earthquake, s, is defined to be the number of top- 
pling events during this earthquake. Then, the force is 
again increased uniformly, until the next site reaches the 
threshold, triggering the next earthquake, and so on. 

Similarly to real earthquakes, which follow the 



Gutenberg-Richter law, the size distribution of earth- 
quakes (which we shall subsequently call "avalanches") 
in this model is found in computer simulations to resem- 
ble a power law n(s) ~ s~ r . The exponent r appears to 
depend on a. However, recent evidence points to a uni- 
versal exponent r ~ 1.8 if only the system size is made 
large enough, at least for a > 0.17 0. Several authors 
argue that there is a critical value of a, around 0.18, be- 
low which the dynamics change qualitatively and become 
dominated by small avalanches [[| |l2| . There is possibly a 
small value of a, around 0.05, below which the power law 
breaks down completely (however, this idea has few 
supporters now), and one recent publication even claims 
that the model is not critical at all for a < 0.25 |l3| . 
This latter claim is based on the finding that the mean 
branching ratio of avalanches is smaller than one, and it 
would be correct if there was only one type of avalanches 
in the system. 

It is the purpose of this paper to clarify some of these 
puzzles and to elucidate the mechanisms that generate 
the observed avalanche size distributions. We will focus 
on the cases a = 0.1 and a = 0.05 and show that the 
avalanche size distribution results from the complex in- 
terplay of several phenomena, including boundary driven 
synchronization and internal desynchronization, limited 
floating-point precision, the slow dynamics within the 
steady state, and the small size of synchronized regions. 
While part of these phenomena have been pointed out 
already in the earlier literature, their combined effect on 
the size distribution of earthquakes and their implica- 
tions for the asymptotic scaling behavior have not been 
clarified so far. It is our prediction that the asymptotic 
avalanche size distribution (in the ideal situation of infi- 
nite floating-point precision) is dominated by avalanches 
of size 1, with the weight of larger avalanches decreasing 
to zero with increasing system size. 

Let us first discuss the case of periodic boundary con- 
ditions, which will be relevant also for the later discussion 
of the system with open boundary conditions. Starting 
with a random initial configuration (most authors choose 
the z values randomly from the interval [0,1)), the dy- 



FIG. 1: Snapshot of the periodic state for L = 100 and a = 
0.1, with an initial width 0.2 of the distribution of z values. 
The force z is visualized by the grey shade, with smaller force 
being darker. 



FIG. 2: Size distribution of avalanches in the periodic state 
for L = 100 and a — 0.1, with an initial width 0.2 of the 
distribution of z values. Dashed line: "Float" precision. Solid 
line: "Double" precision. Cross: Projected result for infinite 
precision. 



namics settle on a periodic attractor, where the sites al- 
ways topple in the same order || ||, |l4| . Often the sites 
are found to topple one by one, without any avalanches 
of size 2 or larger. Periodic states with larger avalanches 
were found for values a > 0.18 0, and they occur also 
for smaller a when the width of the initial distribution of 
z values is made smaller. However, Middleton and Tang 
have shown analytically |jj that a system consisting of 
two sites always settles in a periodic state where none 
of the two sites induces the toppling of the other, lead- 
ing to the expectation that this should also hold for two 
sites that are part of a system with periodic boundary 
conditions in a periodic stationary state. In fact, a closer 
inspection of the periodic states with larger avalanches 
reveals that within an avalanche sites are always lifted ex- 
actly to the threshold value by the force increment they 
receive from the toppling neighbor. At this point, the 
process pointed out by Middleton and Tang comes to a 
halt due to the finite floating-point precision, because 
each site in the system receives during each period the 
same force increment of 4a from its neighbors, allowing 
thus for a periodic state. If the precision was infinitely 
large, the periodic state could only be reached if no site 
induced the toppling of another site. I tested this con- 
clusion by studying a system of linear size L = 100 with 
a = 0.1 and with the initial values of z randomly cho- 
sen from an interval of size 2a. The first avalanche cov- 
ers the entire system, but subsequent avalanches become 
smaller, until a periodic state is reached. Figure |l] shows 
a snapshot of the periodic state, just before the first site 
topples. The pattern visible in this snapshot was gen- 
erated during the first, system-wide avalanche. The site 
that topples first during a cycle has the largest force. The 
bright (but not completely white) sites topple after one 
of their nearest neighbors has toppled. The light grey 
sites (the majority) topple after two of their neighbors 
have toppled, and the remaining ones topple only after 



three or four neighbors have toppled. The force difference 
between two neighbors has therefore peaks at multiples 
of a, as already found by Grassberger for larger values 
of a H . Figure || shows the avalanche size distribution 
in the periodic state, averaged over many different re- 
alisations of the disorder. n(s) is the total number of 
avalanches of size s, divided by the number of sites and 
by the number of cycles. A cycle corresponds to increas- 
ing the force by 1 — 4a. In a periodic state each site 
returns to its original force value at the end of a cycle. 
For larger precision ("double" instead of "float"), there 
are more smaller and less larger avalanches, supporting 
the claim that the floating-point precision limits the de- 
gree of desynchronization. Within an avalanche, all sites 
are lifted exactly to the threshold (within the given preci- 
sion), justifying the prediction that for infinite precision 
all avalanches are of size 1 . 

Middleton and Tang have also shown that if one of 
two sites has a somewhat shorter period than the other 
one (for instance by having a lower threshold value or a 
smaller a), then the two sites settle in a periodic state 
where the site with the longer period always induces top- 
pling of the site with the shorter period. Similarly, sys- 
tems with periodic boundary conditions and frozen dis- 
order in the threshold values or in the values of a can 
have a synchronized state where all sites topple during 
the same avalanche, if the strength of the disorder is not 
too large. (This is shown in M] for disorder in the a val- 
ues. The strength of disorder in the threshold values in 
p5[ was too large to see the synchronization, but I have 
observed synchronization for weaker disorder.) 

When the system has open boundary conditions, sites 
at the boundary receive less force from their neighbors, 
and can therefore be considered as having a smaller pe- 
riod, according to Middleton and Tang B. This explains 
the observed formation of synchronized regions, which 
starts from the boundary and proceeds inward with time, 
apparently according to a power-law la] , which seems to 
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be universal for a > 0.15 After some time (which 

may be extremely long for large L or small a), the sys- 
tem reaches a stationary state with a "patchy" structure. 
We will call a patch a "synchronized region". Within a 
synchronized region, the sites have similar z values and 
topple within a short time of each other, usually through 
a sequence of small avalanches. Regions that are further 
away from the boundary are larger. The regions change 
their shape and size on a time scale that is much longer 
than the timescale of driving, due to large avalanches that 
are triggered close to the boundary. For smaller a, there 
are more and smaller such regions, while for a larger than 
around 0.17, most of the system is dominated by one or 
a few large regions. This paper focusses on the case of 
smaller values of a, where there are several regions of dif- 
ferent sizes, and where the uncertainty about the model's 
behavior is largest. The results reported in the following 
were obtained for a — 0.1 and a — 0.05. Snapshots of 
the stationary state for these two values of a are shown 
in Figure |. This fi gure reveals a striking similarity be- 
tween synchronized regions and a system with periodic 
boundary conditions and a narrow initial distribution of 
z values. Like the system with periodic boundary con- 
ditions, the synchronized regions have different types of 
sites that topple after 1, 2, 3 or 4 of their neighbors have 
toppled. The probability distribution of force differences 
between neighbors has peaks at multiples of a (figure 
not shown, but see for instance Q for larger values of a). 
Just as the pattern of these different types of sites is gen- 
erated by the initial system- wide avalanche in the system 
with periodic boundary conditions, synchronized regions 
are generated by large avalanches, which proceed inward 
from the boundary. Sites that have participated in the 
same avalanche have z values within a limited range, and 
appear therefore as a patch of a given grey shade in Fig. |^. 
They remain in a periodic state where the sites topple in 
the same order for a long time. This can best be visu- 
alized by viewing the system after each cycle, i.e., after 
a force increment of 1 — 4a. The inner part of the sys- 
tem remains unchanged for many cycles. The mean time 
for which a section of the system remains unchanged in- 
creases with increasing system size, with decreasing a, 
and with increasing distance from the boundary. If a 
section remains unchanged after a cycle, each site has 
toppled once and has received a force increment of 4a 
from its neighbors, just as in the system with periodic 
boundary conditions, and often sites are lifted exactly to 
the threshold. The periodic behavior of a synchronized 
region is terminated when an avalanche enters it from 
outside and reshapes the pattern of z values. 

These similarities between synchronized regions and 
systems with periodic boundary conditions (and narrow 
initial width of z values) suggest that the size distribu- 
tion of avalanches is sensitive to the floating-point pre- 
cision. This is illustrated in Fig. |3| In order to demon- 
strate the importance of the distance from the boundary, 




FIG. 3: Snapshots of the stationary state for L — 200 and 
a = 0.1 (top), L = 100 and a = 0.1 (center) and L = 100 
and a = 0.05 (bottom). The force z is visualized by the grey 
shade, with smaller force being darker. 

avalanches triggered within the inner part of the system 
and close to the boundary were monitored separately. In 
addition to using two different types of precision ( "float" 
and "double"), a third simulation was performed, where 
the z value of sites that were lifted exactly to the thresh- 
old by the force increment received from a neighbor was 
slightly decreased (by 10 -15 ) in order to prevent such 
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FIG. 4: Size distribution of avalanches for L = 100 and a = 
0.1., averaged over approx. 200000 cycles. Left: Avalanches 
triggered at least at a distance L/4 from the boundary. Right: 
Avalanches triggered at most at a distance L/8 from the 
boundary. Long dashed line: "Float" precision. Solid line: 
"Double" precision. Short dashed line: Curve obtained if sites 
that are lifted exactly to the threshold do not participate in 
the avalanche. 




FIG. 5: Size distribution of avalanches for L = 200 and a = 
0.1 (left) and L = 100 and a = 0.05 (right), averaged over 
approx. 200000 cycles. Only avalanches triggered at least at a 
distance L/4 from the boundary are considered. Long dashed 
line: "Float" precision. Solid line: "Double" precision. Short 
dashed line: Curve obtained if sites that are lifted exactly to 
the threshold do not participate in the avalanche. 

sites from participating in the avalanche. While for the 
lower precision only 45 percent of all sites in the inner 
part of the system topple in avalanches of size 1, this per- 
centage increases to more than 80 percent with "double" 
precision. In the outer region, there is no such sensitivity 
on the floating-point precision, due to the smaller size of 
the regions and their shorter lifetime. 

With increasing systems size and decreasing a, the sen- 
sitivity to the floating-point precision increases, as shown 
in Fig. [jj. With "float" precision, only 25 percent of the 
sites in the inner region topple in avalanches of size 1 for 
L = 200 and a = 0.1, while this percentage is 75 for "dou- 
ble" precision and 94 when sites that are lifted exactly 
to the threshold do not participate in the avalanche. For 
even larger system size (and the same value a = 0.1), 
on can expect the latter percentage to approach 100, 
implying that virtually all topplings in the inner part 
of the system occur in avalanches of size 1. A similar 



trend is observed for decreasing a. The curves in the 
right part of Fig. || are not smooth due to the extremely 
long time scales for rearrangement of the synchronized re- 
gions. Their cutoff occurs at smaller s than for a = 0.1, 
probably due to the small size of synchronized regions. 

Taking all these observations together, we can conclude 
that the observed avalanche size distributions result from 
the superposition of larger, synchronizing avalanches 
triggered near the boundaries and smaller avalanches 
that occur within the synchronized regions and that 
should turn into single-flip avalanches after a few cycles, 
if only the floating-point precision was infinite. With 
increasing system size and decreasing a, most of the sys- 
tem stays for an increasingly long time in the periodic 
state with single- flip avalanches. The weight of larger 
avalanches would therefore decrease towards zero in the 
thermodynamic limit of infinite system size, if only the 
floating-point precision could be made infinite. The sys- 
tem is thus not critical in the conventional sense that 
it has a scale-invariant avalanche size distribution, with 
large and small avalanches being generated by the same 
mechanism. Instead, the complex interplay of synchro- 
nization, desynchronization, limited floating-point preci- 
sion, and a nontrivial size distribution of synchronized 
regions generates a broad and power-law like avalanche 
size distribution for parameter values typically used in 
simulations. Similar effects might be at the bottom of 
many apparent power laws in nature. 
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